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Abstract —This work presents a statistical analysis of a class 
of jointly optimized beamformer-assisted acoustic echo cancelers 
(AEC) with the beamformer (BF) implemented in the Gener¬ 
alized Sldelobe Canceler (GSC) form and using the least-mean 
square (LMS) algorithm. The analysis considers the possibility 
of independent convergence control for the BF and the AFC. 
The resulting models permit the study of system performance 
under typical handling of double-talk and channel changes. We 
show that the joint optimization of the BF-AFC is equivalent to 
a linearly-constrained minimum variance problem. Hence, the 
derived analytical model can be used to predict the transient 
performance of general adaptive wideband beamformers. We 
study the transient and steady-state behaviors of the residual 
mean echo power for stationary Gaussian inputs. A convergence 
analysis leads to stability bounds for the step-size matrix and 
design guidelines are derived from the analytical models. Monte 
Carlo simulations Illustrate the accuracy of the theoretical models 
and the applicability of the proposed design guidelines. Examples 
include operation under mild degrees of nonstationarity. Finally, 
we show how a high convergence rate can be achieved using a 
quasi-Newton adaptation scheme in which the step-size matrix is 
designed to whiten the combined input vector. 

Index Tenns —Acoustic echo cancellation, adaptive filtering, 
beamforming, generalized sidelobe canceller, statistical analysis 


I. Introduction 

A coustic echoes arise in hands free communications 
when a microphone picks up both the signal radiated 
in a direct path by a loudspeaker and its reflections at the 
borders of a reverberant environment. Acoustic echoes tend 
to degrade intelligibility and listening comfort |[T], . Mod¬ 

ern solutions incorporate adaptive echo cancellers. However, 
typical room reverberation times require adaptive acoustic 
echo cancelers with very long responses m, 0. Also, signal 
contamination by speech from other talkers, noise and their 
reflections in the acoustic environment make it difficult to 
obtain fast convergence and satisfactory echo cancellation with 
such long cancelers m-ni. Moreover, conventional acoustic 
echo cancellation also requires a complex control logic to 
avoid divergence during double-talk periods El, 0. Very 
few studies consider the adaptation during those periods. A 
recent work fSl proposes the use of blind source separation 
techniques. Though promising, such technique still lacks com¬ 
putationally efficient solutions. 
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Assuming it is possible to estimate the direction of arrival 
(DOA) of the desired speaker, spatial filtering (beamforming) 
can help attenuate interfering signals in other directions than 
the desired one. Beamformers (BFs) have limited echo sup¬ 
pression capacity due to limits in the array directivity El and 
the large number of microphones necessary to suppress all 
reflections outside the desired DOA ifTOl . 

Acoustic echo cancellation solutions in which BFs and 
acoustic echo cancelers (AECs) have complementary functions 
have raised a lot of interest recently lfTTII - ll2Tl . BFs and 
AECs contribute by different means to reduce the residual 
echo. Hence, using both techniques in a synergistic way can 
improve the acoustic echo cancellation performance m- 
1221. BFs and AECs are usually combined by means of two 
basic structures mil . l23l . The AEC first structure (AEC- 
BE) employs one AEC per microphone 1191 - 11211 . The BE 
then processes the AEC outputs for spatial Altering. It re¬ 
quires several long AECs, leading to very high computational 
costs HU. Moreover, signals outside the desired DOA must 
be treated as double talk, complicating the design. The BF 
first (BF-AEC) structure does the spatial filtering first, leaving 
basically the echo in the desired DOA to be canceled by a 
single AEC mi-iisi. This structure presents a significantly 
lower computational complexity when compared to the AEC- 
BE structure, even considering that the BF impulse response 
adds to the length of the response to be identified by the 
AEC El]. However, as a single AEC has to cancel echoes 
arriving at many microphones and its desired signal is affected 
by the BE state, the plant identification model is not valid. 
Therefore, previous theoretical work has to be used carefully 
when this structure is studied. In addition, since the AEC 
solution depends on the BF state, an abrupt change in the 
desired DOA can lead to a degraded performance until the 
AEC tracks the new solution. 

Alternative structures that have been proposed include the 
use of polynomial approximations in delay-and-sum beam- 
formers Il24ll . ||25] . the Transfer-Eunction Generalized Side- 
lobe Canceler (TE-GSC) Il26l - ll29l . AEC sub-modeling 
mutually exclusive adaptation of the BE and AEC m and 
wave-domain Altering ll^ . ll^ . 

Optimization of BE-assisted acoustic echo cancellation sys¬ 
tems can be based on different performance surfaces, depend¬ 
ing on how the BF and the AEC are optimized. One may 
define the beamformer performance surface from its own local 
error E2 or use a joint optimization scheme 01, El in 
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which the global cancellation error is used to jointly optimize 
the BF and the AEC. The joint optimization scheme was 
first proposed in ifT^ . It was later applied to a robot speech 
recognition system Qa. Joint BF-AEC optimization leads to 
an optimal solution with better echo cancellation performance 
than separate BE and AEC optimizations fT2\ . 

Despite the possibilities of combined BE and AEC acoustic 
echo cancellation systems, we find only few analyses of their 
transient behavior in the literature. The AEC-BE structure 
has been studied for the acoustic echo cancellation problem 
in mi-iEi] and for the acoustic feedback cancellation in 
ll34l . A stochastic model has been derived using the power 
transfer function method for the case of a fixed BE, where 
just the AEC is adapted. More recently, the transient behavior 
of a system where a direct-form BE and an AEC are jointly 
adapted using equal and fixed step-sizes was analyzed in El, 
M- The derived analytical model was shown to accurately 
predict the adaptive system behavior and corroborated previous 
experimental findings that the same cancellation performance 
of a single-microphone AEC can be achieved with a shorter 
AEC when the possibility of spatial filtering is available llT5l . 
The model, based on the equivalence to a coventional Linearly 
Constrained Minimum Variance (LCMV) optimization, allows 
the use of previous analytical results M, El- 

Adaptive LCMV beamforming may be implemented in 
many different forms and by using different algorithms 1^ . 
El-ESl- The direct and GSC forms are equivalent in that 
both lead to the same optimal solution mu. Eor some algo¬ 
rithms and under specific conditions they are equivalent even 
in their transient behavior 1^ . 1451 . Both forms 

tend to have comparable computational complexities for a 
small number of constraints. Flowever, the GSC form offers 
greater design flexibility due to the possibility of choosing the 
block matrix. Good choices may lead reduced computational 
complexity 1^ p. 31]. Also, robust GSC implementations 
with an adaptive bock matrix have been proposed to account 
for small changes in the desired signal DOA ITSl . I46l - l49l . 
Therefore, it is of interest to study the behavior of the GSC 
form of the BE-AEC structure. 

This work extends the analysis in IITtI . ITSl to the study 
of the transient behavior of the jointly optimized BF-AEC 
structure in the GSC form. We formulate the joint optimization 
as a single constrained optimization problem, what simplifies 
the statistical analysis. Moreover, the analysis incorporates 
the case of a positive-definite step-size matrix ll50l - ll54l . The 
incorporation of this extra flexibility to the model is particu¬ 
larly interesting for BE-assisted echo cancelers, as their AEC 
adaptation control logic stops AEC adaptation during double- 
talk periods la, Q, while the BE continues adapting using 
a Reference Signal Based (RSB) structure 1551 with the AEC 
output as the reference signal. The problem of designing an 
adaptive filter with step-size matrices was studied in Il50l - ll54ll . 
An exponential model for the echo channel and information on 
the room reverberation time were exploited in 15^ to design 
a step-size optimized algorithm. In Il54ll . it was shown that 
LMS algorithm with a step-size matrix is equivalent to the 
classical LMS algorithm in a transformed space. The same 
idea is used in our convergence analysis. The analytical model 


derived in this paper allows the study of the echo canceler 
behavior including echo-only periods, when AEC adaptation 
is slower, double-talk periods when only the BE is adapted, 
and periods after channel changes when fast AEC adaptation 
is required ||6l, Q. 

The main contributions of this paper are: 

(i) The formulation of the jointly optimized BF-AEC im¬ 
plemented in the GSC form as an LCMV-based GSC. 
This signal model can be used to design the conventional 
LCMV-based GSC without loss of generality. Previous 
theoretical results show that the behavior of the GSC 
can be studied from the direct form when adaptation 
uses a single step-size, feasible quiescent solutions and 
blocking matrices have orthonormal columns [381, ED- 
Hence the analysis can also be used to design the BE- 
AEC and conventional LCMV implemented in the direct 
form using a scalar step-size generalizing the analysis 
in El; 

(ii) Incorporation of a step-size matrix. AEC adaptation 
control logic demands the adaptation of the AEC and BF 
with different step-sizes during different adaptation sce¬ 
narios (double-talk, channel changes, tracking, etc) Q, 
Q. Hence, a novel analysis capable of predicting the 
transient behavior during different control logic states 
(different step-sizes) is of undisputable practical rele¬ 
vance. The analysis model uses a positive-definite step- 
size matrix 

Using the proposed formulation, we derive a statistical model 
of the behavior of the BF-AEC system implemented in the 
GSC form with a positive-definite matrix step-size. The model 
also allows the derivation of a high convergence rate algorithm 
based on a quasi-Newton adaptation scheme in which the step- 
size matrices are designed to whiten the combined input vector 
to accelerate convergence. 

This paper is organized as follows. Section II formulates the 
problem addressed. Up to Section II-C the material is basically 
the same as in ifTSl and is necessary to establish the notation 
used in the rest of the paper. Section II-D introduces the GSC 
formulation for the problem studied. Section III describes the 
analysis structure that allows the analysis of the adaptation 
using different step-sizes and the quasi-Newton algorithm 
using the same mathematical framework. Section IV describes 
the adaptive solution. Section V derives the statistical model 
for the adaptive solution. The statistical model convergence is 
analyzed in Section VI. Based on the results in section VI, 
the new quasi-Newton adaptation is derived in Section VII. 
Section VIII validates the proposed model using simulation 
examples. Einally, conclusions are presented in Section IX. 
In this paper, plain lowercase or uppercase letters denote 
scalars, lowercase boldface letters denote column vectors and 
uppercase boldface letters denote matrices. 

II. Problem Eormulation 

Fig. m shows the BF-AEC structure with M echo impulse 
response vectors hm of length Ah, M microphone signals 
XTn\n], one adaptive wideband beamformer composed of M 
filters bm[n] of length A^bf and an adaptive AEC filter h[n] 
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of length Naec- We assume impulse responses constant 
and stationary signals for mathematical tractability HI pp. 
348-351]. The analysis for a time variant echo path becomes 
specially challenging in this case even for the simple random 
walk system nonstationarity model Q, ||5l. This is because a 
time variant loudspeaker-enclosure-microphone (LEM) model 
would lead to a nonstationary beamformer input signal. More¬ 
over, the statistically independent increments to the channel 
response vectors hm due to the random walk model would 
be time-correlated by the BE filters. This would render the 
analysis too complex even for such simple nonstationarity 
model, making it very hard to study fundamental properties 
of the algorithm behavior. The study for nonstationary input 
signals requires a specific model for the input nonstationarity. 
To the best of our knowledge, there is no generally accepted 
model for signal nonstationarity. On the other hand, model 
predictions derived under stationarity assumptions can still 
show tendencies of the algorithm behavior for reasonably 
small degrees of nonstationarity Up. 595]. Simulation results 
in Section IVlll-CI will illustrate that this is the case for the 
present study. 

It has been conjectured that the spatial filtering realized by 
the BE reduces the required AEC length, as compared to the 
conventional finite impulse response (EIR) AEC structure llTSl . 
Hence, our analysis considers the possibility of an AEC shorter 
than the LEM impulse responses by admitting Aaec < ^h- 



Fig. 1. BF-AEC system configuration in the direct-form structure (22). 


A. The Beamformer Input Vector 

Each of the M LEM impulse responses hm, m = 
0,..., M — 1, models the transmission of the far-end signal 
u[n] from the speaker to one of the M microphones. The 
adaptive wideband beamformer is composed by M EIR filters 
with impulse responses bm[n], m — 0,..., M — 1, each of 
length TVbf ll57l . The echo signal at the mth microphone is 
given by H] 

em[n] = h^Ui\n] (1) 

where 

tth[n] = [u[n\,u[n - 1], • ■ • ,u[n- (iVh - (2) 


Grouping the LEM responses as columns of the matrix 

if = [/lo hi ■■■ Hm-i] (3) 

and defining the echo snapshot vector as 

es[n] = [eo[n],ei[n],-- - ,eM-i[n]f (4) 

([T]i leads to the linear mapping 

es[n] = iT^Mh[n]. (5) 


The mth microphone signal Xm [n] is the sum of a near-end 
signal rm[n] and an echo emln]- 

Xm[n\ = em[n] + rm[n\, m = 0,...,M-l. (6) 


Each signal rm[n] is composed of local speech, local inter¬ 
ferences and random noise. We define the microphone array 
snapshot [n] as the vector composed by all Xm [n]: 


x^[n\= a;oN,a;iN, • • • ,XM-i[n] 

Then, combining (01, (|6]) and (|7]i yields 


n T 


(7) 


a:s[n] = es[n] -f rs[n] 

where rs[n] = [ro[n],ri[n], • • • ,rM-i['n]]^ is the near-end 
signal component snapshot. 


We now define the extended far-end sample vector as 
u[n] = [u[n],u[n - 1], • • • ,u[n - (Ah -f Agp - 2)]]^ (8) 


where the dimension of u[n] is the length of the convolution 
of hm[n] and bm[n]. Then to express the microphone array 
input signals (the echo signals) as functions of u[n] we rewrite 
(|5]l as 


eJn 


-k] = 


JMxk 




0 


M X Abf~ (fc+ 1 ) 


u[n] (9) 


where 0 i\[-ixN 2 denotes the null matrix with dimension Ai 
lines and N 2 columns. Then, defining the M.Nsp x 1 stacked 
echo vector 


eN = [ej[n],ej’[n-l],--- , - (Abf - 1)]]^ (10) 


we can write 


e[n] = 'HT' u[n] 


( 11 ) 


where 



1 OlxM 1 ■ ■ 

1 ®A^bf —lxA(f 

'H = 

H 1 H 1 

1 H 


^Abf —IxAF ^Abf — 2xAf 

1 


is the Ah + Abf — lx M.Nbf modified echo channel matrix. 
Note that e[n] contains the echo signals for the time window 
corresponding to the length of the BE impulse response. 


Using (|6ll, (I 2 I), ([Tol l and (fTTT l. and defining the M.Abf x 1 
near-end vector component (without echo) as 

T-b[n] = [rf [n],rf [n - I],-- - ,rj'[n - (Abf - 1)]]^ (13) 

and the M.Npp x 1 combined beamformer input regressor 


is the LEM plant input vector. 
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Xi,[n] = xj'[n],xj'[n - 1], ■ 


[n — (-/Vbf “ 1)] 


n T 


(14) 


we write the beamformer input vector as 


Xh[n] = n u[n] + rh\n]. 


(15) 


B. The Residual Echo 

Define the vector \n] of the £th components of all vectors 
bm[n], m = 0,..., M — 1, at time n as 


K [n] = 




= 0,..., A^bf - 1. 


We then write the beamformer output y[n] as 

A^bf—1 

2/W = xj[n-i]b,,[n\. 


e=o 


(16) 


Now, defining the stacked beamformer weight vector 

„ „ iT 


b[n] = 




Ms 


we can write y[n] as the inner product 
y[n\ = x^[n\b[n\. 

Next, defining the AEC weight vector 


h[n] = 


hoN.^iN,-- 


) ^tVAEC — 



(17) 

(18) 

(19) 


and the AEC input vector 

= [u[n],u[n - 1], • • • ,u[n- {Naec - 1)]]^ (20) 

we can write the AEC output as 

y[n\ = [n\uf^[n\. (21) 

Using (fTSl l and (|2T]) we write the residual echo d[n] as the 
inner product 

d[n\ = —uj [n]h[n\ + a;|f [n]b[n]. (22) 


III. The Analysis Structure 

With the problem formulation presented in Section [H] we 
can define an analysis problem that corresponds to the study of 
a single GSC structure that combines the beamformer and the 
AEC adaptations. To this end, we define the Waec+-^-^bfX 1 
stacked input vector 

s[n]= [-U(^'^[n\,x^'^[n]]'^ (23) 


and, from (fTTI l and (fTOl l, the stacked coefficient vector 


iu[n] = 


h 


(24) 


Then, we can write the residual echo d[n] as the inner 
product 


Interestingly, input vectors u[n] in ([8]) and u^[n] in ( l20l i are 
related by 

“hW = [ -^YaeC OiVAEcXiJVh + JVBF-JVAEC-l) ]“N (26) 

where we use the notation Ik io denote the K x K identity 
matrix. Hence, ([I5]l and (|26 T i permit to write s[n] in (l23t as a 
function of the input vectors u\n] and r\,[n]. Equation (l26l l also 
allows to study the algorithm performance for A^h > ^aeCi and 
thus verifies the possibility of reducing Waec by increasing the 
number of microphones. 


A. Performance Surface 

The mean output power (MOP) performance surface J is 
defined as the mean value of d‘^[n] conditioned on iu[n] = w. 
Erom (|25 T i. 

J = = w} = E 

= w^RssW. (27) 


where R^s = i?{s[n]s^[n]} is the input autocorrelation ma¬ 
trix. A set of Nf linear constraints on the beamformer coeffi¬ 
cients implements the spatial filtering. Usually, an MNbfxN/ 
constraint matrix C and an Nf x 1 response vector / jointly 
define the frequency response in the desired DOA 1^ . 1^ . 

To formulate the linear constraints as a function of the 
combined coefficient vector, we define the extended constraint 
matrix Cl 


Ce 


X Naec 



(28) 


Finally, the joint optimization problem can be formulated as 


lUopt = arg min 


w^RssW 


subject to Cg w = f 


(29a) 

(29b) 


and the optimal solution is given by ll^ tUopt = 

r-Me [c^R-M,y'f. 


B. Implementation using the GSC Form 

In the GSC form 1^ . the dashed square of Fig.[T]is replaced 
by the dashed square of Fig. H Feasible solutions to (|2^ are 



Eig. 2. BE-AEC system in the GSC configuration. 


decomposed as ll38l 


d[n\ = s^[n]w[n\. 


(25) 


w = q^- Be-ip 


(30) 


This simple model will permit to relate the study of the 
BE-AEC structure to that of the LCMV problem. 


where is any feasible solution to (I29bl l. is a full column- 
rank {Naec + MN^f) x -dimensional blocking matrix 
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orthogonal to Ce (CjBg = 0), ip is an -dimensional 
vector and = Naec + MNsp — Nf. The minimum norm 
solution to (I29bl i is 

q, = Ce(CfCe)-V- (31) 

1) Optimal Solution: As CjBe = 0, m in dSOl l satis¬ 
fies (I29bl l for any r/;, and (|2^ becomes an unconstrained 
optimization problem in xp with solution Bdl 

Ip t = a-Tguiinql- 2ip'^BJR^^qe + ip'^Rhiocip 
V 

(32) 

where i?bioc = -B^-Rss-Be denotes the blocked input autocor¬ 
relation matrix, and from ( l30l l 

«topt = Qe - Belpepf (33) 

Defining the cost function of ( |32] | 

C{ip) = qjRssqe - ‘2tp^BjRs,qe + ip'^Rhioc'ip (34) 
its gradient with respect to “ip is 

V^C(-0) = -2BjR,,{q^ - Be-iP) (35) 

Setting (l35l l equal to the null vector yields HTI 

^opt = Rbioc^J R.,qe- (36) 


IV. The Weight Adaptation Equation 
To obtain a model flexible enough to allow the study 
of the system performance with independent BE and AEC 
adaptations we choose the following block diagonal form for 
the blocking matrix 


Be = 


-^Naec 

0 M Abf X Naec 


Unaec X (MNbp — N^) 

B 


and split 




(37) 


(38) 


where [ipf^]i = [ip]i, i = .,Naec and [ip,,]i = [ip]i+NAEC’ 
i = 1,., MNaec — Nf. The same block matrix structure has 
been used in 1(141 . ifTSl for the implementation of the GSC- 
based BE-AEC acoustic echo canceler. 

Using (l30l l in (l35l l and noting from (l25l l that R^sW = 
i?{s[n]d[n]} we have 

V^C{iP) = -2BjE{s[n]d[n]}. (39) 

Splitting the gradient vector according to (l38l l yields 

V^Ci^)= [vl.C{xP),VlC{xP)\^ 

where, from (l23l l, (iJTl l and 


= 2E{uf^[n]d[n]} (40a) 

\7^^C{ip) = -2B^ E{x^,[n]d[n\}. (40b) 

Additionally, setting q^ = [OixAaec 9^]^ and us¬ 
ing dJTl l and OTl l in (l30l l yields 


w = 




(41) 


where q = C{C^C)-^f. 

Comparing (l24l i and (l4Tl i we conclude that tpf = h and 
b = q Bipii. Hence, the steepest-descent algorithms for h[n] 
and '0b[n] with the gradients in (I40al and (I40bl respectively 
are 

h[n -I- 1] = h[n] - paecE{ u{^[n]d[n]} (42a) 

ipi,[n -I- 1] = 'ipbln] + pbfB"^E{ xb[n]d[n]} (42b) 

where /taec and pbf are the step-size parameters. Note 
that (l42l i is different from the steepest descent algorithm for 
■0[n] unless /taec = Mbf- This also makes this analysis 
different from ESI by using the equivalence derived on ll38l . 
ED. However, this extra degree of flexibility is necessary to 
analyze the behavior of the BE-AEC system under different 
control logic states that usually act on (/taeCj/^bf) to avoid 
divergence. The stochastic approximations of (I42ab and (I42bb 
yield 

h[n + 1] Fi h[n] - pAEGUb[n]d[n\ (43a) 

^Pl,[n -I- 1] ~ 'iPbM + p.BFB'^Xb[n]d[n]. (43b) 

Implementation of (l43l l has almost the same computational 
cost of the separate implementation of an LMS implementation 
of an AEC and a BE demanding only an extra subtraction 
in the computation of d[n]. It also requires only one extra 
memory allocation to account for the second scalar step-size. 
Despite its simplicity, EUl can model the BE-AEC system 
behavior under most control logic states. Implementation of 
(l43l l is shown in Fig. |3] 



Fig. 3. Implementation of g3) 


Finally, the recursive weight update equation is obtained 
defining the diagonal step-size matrix 


M.= 


PaecIna, 


0 


Naec X (MN^F — Nf) 


N^p — N f) X Naec l^EF^^ M Nef — N f 

then (HD can be written as 

■0[n -I- 1] = ip[n] + AiBjs[n]d[n]. 


(44) 


(45) 


Note that (l4Tl i has the exact same behavior of (l45l l. which can 
be used to study the performance of the practical implemen¬ 
tation. 

In the following we perform the analysis of an even more 
general form of (l45l l. in which the only requirements on AT 
and Be are that AT is symmetric positive-definite and Be 
is a full column-rank matrix that satisfies CjBe — 0. The 
typical implementation described above will correspond to a 
particular case of the more general analysis. 
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A. Weight Error Vector 

Define the weight error vector v[n] = w[n] — ritopt- 
From d^ . 

■wW = ge - - (q, - -BetAopt) = (46) 

where 

-d[n] =ip[n]- (47) 

denotes the weight error vector of the unconstrained hlter 
conditioned on and q^. From (|46] |. v[n\ is in the range of 
JBe- Hence, v[n] is completely determined by 'i9[n] conditioned 
on Be- We then study the behavior of 'i9[n]. 

Subtracting from both sides of (l45l l. using (l46l l with 
(l3?t and (l30l l we obtain a recursive update equation for 

i9[n + l] =ilN^ - MBJs[n\s^[n]B^)-d[n\ 

rp rp ('^O) 

+ M.B^ s[n]s [nJiUopt- 

V. Statistical Analysis 

A. Simplifying Assumptions 

We now study the behavior of BF-assisted GSC-form echo 
canceler using (l45l l under the following typical simplifying 
assumptions required for mathematical tractability B 

A1 s[n] is a zero-mean Gaussian vector; 

A2 u[n] and r[n] are statistically independent; 

A3 i?ss is positive-definite and both Ce and B^ have full 
column rank; 

A4 The statistical dependence between s[n]s^[n] and 
■0[n] can be neglected; 

A5 The DOA does not change during adaptation. 
Though not always valid in practice, these assumptions make 
analysis viable and frequently lead to results that retain suf¬ 
ficient information to serve as reliable design guidelines IS 
p. 315], Ida, m. Simulation results will confirm their 
reasonability for this analysis. A1 simplifies the evaluation of 
fourth order moments of s[n]. These moments are dependent 
on the distribution of s[n], and the Gaussian distribution 
combines the advantages of being a good model for several 
physical processes and simplifying the required mathematical 
derivations. A2 is physically reasonable, as u\n] and r[n] are 
generated at different sides of the communications channel by 
independent speakers. A3 is reasonable in practice, as s[n] 
always has some uncorrelated noise component and both Ce 
and Be are under reasonable control of the designer. A4 is 
required to estimate moments involving the input signal and 
the weight vector, as the statistical distribution of the latter is 
unknown. This assumption is in fact less restrictive than the 
usually employed independence assumption, which requires 
s[n] and tp[n] to be independent, as discussed in detail in ll58l . 
A5 is employed for mathematical tractability and because the 
main goal of the present analysis is to determine fundamental 
properties of the adaptive system. 

B. Mean Weight Error Vector Behavior 

Taking the expected value of (l48T l under A4 and using (l30] l 
and (l36l l leads to 


E{B[n + 1]} = {In^ - MRuoc)E{B[n]} (49) 

since 

B{Bg s[n]s [njlrUopt = -Bg -Bss^e ~ -^bloc-Bbloc-®e -Bss^e 
= Ojv„xi- (50) 

Hence, the mean weights converge asymptotically to the 
optimal solution if all eigenvalues of In^ — ATBbioc are 
inside the unit circle. In this case, (l45T l results in asymptotically 
unbiased solutions in the mean. 

C. Mean Output Power (MOP) 

To determine the MOP we use (l25l l with w[n\ = t;[n]-|-iUopt 
and (l46l l. Defining B^^[n] = E{'d[n]'d"^[n]} we obtain 

J[n] = E{{wopt - Bei9[n])^s[n]s'^[n](i(Jopt - Bei9[n])} 

— '-^min T tr(B^^ [ri] Bbloc) (51) 

where we have used (fSOl l and A4 to obtain the second line. 
A recursive expression for the x matrix i?ij^[n] is 
derived in the next section to complete the model (ISTT i. 

D. Correlation Matrix of 

Post-multiplying (l48T l by its transpose, taking the expected 
value, using A1-A5 and (fSOl l yields 

4-1] — A^.Bbioc-B^^[7z] R,'0'ff^Ti\R.\yiQcJVt 

+ ME^B^ s[n]s^ [n]B^'d[n]'d^ [n]B^ s[n]s^[n]Be| A4 

4- JminTWBblocTW. (52) 

Using Al, A4 and the Gaussian moment factoring theorem, 
the expectation in (l52l i is given by 

E{B^ s[n]s'^ [n]B^'d[n]'d^ [n]Bj s[n]s^[n]Be} 

— 2F^bloc-B^T9 [tt]-Bbloc 4“ .Bbloc^(-Bbloc-B^T9 [^]) ■ (53) 

Finally, substituting (l53T l into (l52l i yields 
R''d'd\jt 4“ 1] —[n] A^RyhlocR''&'&\jt\ [77.]Fibloc«^^ 

4“ [<7min 4“ tr(F?bioc-Bi9Tt [rt])]ATF^bioc-^^ 

4“ 2A^.Bbioc-B^^ [n] Fibioc• 

(54) 

Equation (l54l i completes the MOP model in (ISTT i. 

VI. Convergence Analysis 
Classical convergence analysis of dSli would project 
R§s [n] into the eigenspace of Bbioc and study the convergence 
of the diagonal entries of the transformed matrix El. The 
presence of A4, however, requires a different approach. As 
ATBbioc ^ -RbiocA4, (l54l l is not entirely diagonalizable by the 
same projection ||59] p. 558]. Nevertheless, it is still possible to 
diagonalize both A4 and Bbioc through contragradient diago- 
nalization ll59] p. 4651. ll60l p. 466]. As A4 is positive definite, 
Cholesky decomposition yields AT = LL^ with L non¬ 
singular. Then, we can transform the vector space into ^[n] = 
L-^B[n], B^^[n] = B{|[n]|'^[n]} = L-^Rm[n]L-'^ and 

B^d = B{T^Bjs[n](i^Bjs[n])^} = L'^RuocL. (55) 
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Hence, pre-multiplying ( l54l i by i post-multiplying by i ^ 
and using tT{R^^[n]R^od) = ii{Rm[n]Rb\oc) yields 

[tt “b 1] — [tt]-^mod 

“b -^mod['^min “b tr( [tz] ^^niod)] “b 2^];jiQd [tz] ^^ixiod (56) 


where i?mod is symmetric and positive definite. Hence, it is 
diagonalizable as iimod = QAQ^ with Q^Q = 1^^ and 

A = diag(Ai, A 2 ,..., Aw.^). (57) 

Pre-multiplying (l5^ by and post-multiplying by Q yields 


-p 1] = 0^55[n] - Aa^ 44 [n] - [n]A 
+A(J^in + tr(lH 45 [n]A)) + 2A?H^^ [n]A (58) 

where 9i^{[n] = Q'^R^^[n]Q. 

is an autocorrelation matrix. Then < 

[£Hcc[n]],,, > 0 EU p. 251], EH,’ and 
convergence of (l58l l can be studied observing only the diagonal 
elements of Let i>'[n] be the vector of diagonal entries 

of and A = [Ai, A 2 ,..., Aat^]^ be the vector of the 

eigenvalues of i?mod- Then, from (|58] | 


[v>[n + l]]i 
and 


(1 ~ Ai)^ -b Xi + Ai(A^t'[n] -b Jmin) 

(59) 


v>[n -b 1] = ^u[n] -b JminA (60) 

where A^i/[n] = tr(A9^55[n]), 

$ = diag(pi,p2,---,PJV,^) + AA^ (61) 

and pfc = (1 - Afc)^ -b A|. 

The matrix $ is symmetric and positive definite, as for any 
nonzero vector a we have 

= a^diag(pi, p 2 , ■ ■ •, PN^)a + a^XX^a 

= ^ ((1 - Xkf + xl) ([ajfe)^ + (A^a)^ > 0. 

fc=i 

The solution to EQI* is 

n—1 

iz[n] = + J„,i„ ^ #^ A. (62) 

f=o 

Using E3 we now study the stability conditions and the 
steady-state behavior of (l45T l. 


A. Mean Weight Error Revisited 

Pre-multiplying ( |49] | by L~^ the mean weight error recur¬ 
sion is transformed to ED 


1,..., N^, of $ EH- From Gershgorin’s theorem IMl , 

Xif>^ < Pfe + Afc 

i=i 

< 1 — 2Afe "b 2A^ -b Afctr(i?n,Q(j), Vfc (64) 

and ED is stable if A.^^ < 1 for all k. Then, (l64l i leads to the 
sufficient condition 

~ 2Afe "b 2A^ "b Afctr(i?u,Q(]) < 0, Vfc (65) 

which implies that Xk ^ 0 and 

2 max{Afc} -b tr(iiniod) < 2. (66) 

In most practical cases, a reliable estimate of the eigenvalues 
of Rmod is not available a priori and the upper bound in (l66l l 
can not be used. However, using the inequality 

max^Afe} ^ tr(J^j];iod) 

it is possible to derive a tighter upper bound E5l , Il66l 

2 

tr(i?mod) = X! ^ Q ■ 

k^l 

In the particularly important implementation using (EDl, 
Be is given by ED and AT by (l44l i. Hence, if we write 
tr(iiniod) = tr(L^i?bioc^i) = tr(A4i?bioc) with the matrices 
in the partitioned form using (l2Tt . ETl l. (l44l l and ED yields 

tr(i?mod) = AtAECtr(iluj^uJ + /rBFtr(B^ilxbXbB) (68) 

where R^^^^ = E{uj^[n]uT[n]} and R„^^^ = E{xb[n]x'^[n]}. 
Hence, ( l67b becomes 

2 

/iAEctr(Buf,uj) + pBFtr(B^iixbXbB) < -. (69) 


C. Excess MOP 

Using tr(i2bibiBbioc) = ^(^^^[njBniod) = tr((H 55 [n]A) = 
X^u[n] in (fSTI) we write the MOP as a function of A and u[n\ 

J[n] = Jmin + X^ u[n]. (70) 

Thus, the excess MOP is given by El P. 302] 

Jex[n] = A^iz[n] (71) 


D. Steady-State Excess MOP 

When ED holds, lim„^oo zz[n -|- 1] = iy[n] = i>'[oo] and 
from ED and ED we have 


E{^[n + 1]} = (Ja^ - Bmod)B{^[n]}. (63) 

Hence, analysis can be restricted to the eigenspace of Bmod 
defined in (fSD . 

B. Stability Conditions 

Recursion ED is a state-space equation whose stability 
is determined exclusively by the eigenvalues A^;,, fc = 


which solved for [zz[oo]]i yields 

[tz[oo]]i = (Jex[tX)] “b Jmin) 


1 


(73) 


2 - 2A, 

Using this result in ED as n —00 and solving for Jex[oo] 
yields 

1 Xi 

2 1-Ai 


Jex [cxd] — Ju 


1 1 A 

2 Z^2=l 1- 


(74) 










8 


From (|55]i, fibioc = and A3, i?mod is symmetric 

and positive definite. Hence, its largest eigenvalue is related 
to its largest singular value through max{A} = y/ max{cr} 
where cr denotes the vector of singular values of i?n,od- The 
largest singular value of a matrix is equal to its 2-norm WT\ 
pg.78]. Then, using the Cauchy-Schwarz inequality ll^ pg. 
291] 

max{A} < (^||T'^|| 2 ||i?bioc|| 2 ||i^|| 2 ) "" 

< (||M||2||i^bloc||2)^ (75) 


where both A4 and i?bioc are symmetric positive definite 
matrices. Hence, for max{A^} max{A/{j^,^^} ^ 1, where 
and are vectors containing the eigenvalues of AT 

and i?bioc, respectively, we conclude that max{A} ^ 1 and 
(l74l i reduces to 


t7ex[tX)] 


2tr( J^niod) 

1 ^^^("^niod) 


(76) 


For implementations using (l43T l. substituting dhST l in (l76l l 
yields 


/iAECtr(-Ru|'uJ + MBFtr(.B^.RxbXb-B) 

2 — MAECtr(-RubUb) — MBFtr(S^ilxbXb-B) 

(77) 


Further assuming /rAECtr(iiu(Uj) + AiBFtr(B^i?XbXb-B) < 2 
we have 


Tex [c 


fiAECtl'(-Ru(,uJ + MBFtr(S^ilxbXb-B) 


(78) 


VII. A New Joint Adaptation Algorithm 

The analysis results derived in sections [V] and |Vl] are valid 
for the general weight update equation ( l45T l. At the same time, 
(|45T i in its general form where AT and B^ satisfy only the 
criteria of AT being positive-definite and CjBg = 0 can 
be considered a new adaptive algorithm that allows weight 
updating in directions that do not correspond to the stochastic 
gradient. 

Next, we discuss one possibility of taking advantage of 
the more flexible structure, namely, designing for a faster 
convergence speed 

The simplest way to guarantee asymptotic convergence to 
is set AT = plIn^ in (|45]l where ^ = ^aec = fiBF, 
which is the standard LMS update. However, it is known that 
LMS presents a low rate of convergence when the gradient of 
the performance surface has a low magnitude in the direction 
of at least one eigenvector of iJmod- To alleviate this issue, one 
may use the step-matrix in (l44li or a stochastic approximation 
of the Newton method. The idea underlying quasi-Newton 
methods is to use an approximation to the inverse Hessian. The 
form of the approximation varies among different methods - 
ranging from the simplest, where it remains fixed throughout 
the iterative process, to the more advanced where improved 
approximations are built up on the basis of information gath¬ 
ered during the descent process ll68l . 


A. High Convergence Rate Block Matrix and Step-Matrix Pair 

Theoretical results show that the rate of convergence of 
(l45l l is increased with the reduction eigenvalue spread of 
L^B^R,,B,L a, El, m, reaching its maximum when all 
eigenvalues are equal. This is because the MOP of an adaptive 
filter trained with an algorithm of the LMS family decreases 
over time as a sum of exponentials whose time constants are 
inversely proportional to the eigenvalues of the autocorrelation 
matrix of the filter inputs ifTsl . Hence, small eigenvalues create 
slow convergence modes while large eigenvalues limit the 
maximum step-size that can be chosen without encountering 
stability problems as observed in (HU Go). When all R^od 
eigenvalues are equal, we have 

Q^RmodQ = (79) 

where A = \i{R^od)/N^ ~ 

multiplying ( |79] | by Q, post- multiplying by and noting 
that QQ^ — we conclude that jiimod = In^ - Observing 
.Rmod structure from (fSSl) yields (ji?bioc-T) = I and 
observing that both and -^iibiocT are square matrices, we 
conclude that 

= Qilbioci) = AL-ifibloV ( 80 ) 

Finally, pre-multiplying dSOl) by L and substituting AT = 

LLF yields 

^ ( 81 ) 

J[ooJ 

Direct use of dSTT i would require prior knowledge of second- 
order statistics of s[n] and Jmin- Nevertheless, the transient 
behavior when using (IHTT l is a useful measure of the upper 
bound on the convergence speed. A compromise solution 
would be the estimation of AT every few iterations and the 
use of update (l45l l. In this case, (ISTI) becomes a quasi-Newton 
adaptive filter 11711 . 

VIII. Simulation Examples 

This section presents simulation and design examples to 
verify the accuracy of the derived model and to illustrate its 
use in design. In all simulations, except when explicitly stated, 
the far-end signal was drawn from an autoregressive process 
ARl(ai) given by u[n] = — aiM[n—l]-|-z[n], with 2:[n] a white 
Gaussian noise with variance al such that cr^ = 1. For the 
accuracy tests, adaptation of the coefficients is assumed to be 
done during a single-talk period. Simulations under different 
control logic states are considered in subsection IVIII-FI where 
only the BF is adapted during a double-talk period and the 
convergence of both filters is accelerated after a LEM plant 
change is detected. The LEM plants are designed following the 
procedure outlined in ESI for a uniform linear microphone 
array oversampled by a factor of F = which yields 
spatially correlated exponential impulse responses. 

is the ratio between the temporal oversampling factor and the spatial 
oversampling factor used to generate spatially correlated LEM impulse 
responses (HI Appendix B]. 
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A. Model Verification 1 

The accuracy of the derived model has been verified through 
Monte-Carlo simulations using several different parameter 
sets. To conserve space. Fig. |4] shows a few of these sim¬ 
ulations for M = 2 and LEM responses hi and fi ,2 with 
A^h = 128. The beamformer filters had Nbf = 16 and 
linear phase in the look direction. The AEC length was 
IVaec = 128. The noise variance at each microphone was 
10“^. The theoretical predictions (smooth red curves) are 
in very good agreement with the Monte Carlo simulations 
(300 runs). Values of J[oo] using (l76l l are shown by the red 
horizontal dotted lines. 




2/3 2/30 


[3.9840 X 10“"^, 0.0028]. Eig. |5] shows excellent agreement 
between theory and predictions in both cases. Counterintu¬ 
itively, the results show that a larger convergence speed does 
not necessarily imply a higher steady-state error. 



Fig. 5. Monte-Carlo simulation results for 20 runs with ARl(-0.9) input 
{M = 2, = 500, A^bf = 16 A^aec = — 1)- 





(e) white noise /(/jaeCi/^bf) = (0 white noise /(/^AEC, Mbf) = 

2/3 2/30 


Fig. 4. Proposed model and Monte-Carlo simulation results based on 300 
runs for different far-end signal statistics (M = 2, = 128, F = 4, 

Nbf = 16 Waec = 128) 


B. Model Verification 2 

Consider a unit power first order autorregressive AR1(- 
0.9) far-end signal, 2 microphones, hg and hi with 500 
taps each, generated according to the model in im. The 
desired DOA was assumed orthogonal to the microphone 
array. We assumed the absence of double-talk, and noises 
rp [n] and ri \n] were zero-mean white Gaussian with variance 
10“^. The adaptive BE was designed with A^bf = 16, linear 
phase, and all-pass frequency response with Nf — 16. The 
AEC used A^aec = A^h + A^bf — 1- Fig- 13 shows the pre¬ 
dicted and simulated transient MOP. We tested 2 scenarios: 
[/iaeCjFbf] = [2.6191 X 10“^, 0.0262] and [/iaeCjFbf] = 


C. Model Verification 3 

The model in this paper is derived under stationarity as¬ 
sumptions for the input signals. Nevertheless, for nonstationary 
signals, it preserves sufficient information about the adaptation 
process to derive useful design guidelines. It is important 
to stress the fundamental difference between design guide¬ 
lines and design rules. Design guidelines are not rules that 
should be followed to design the system with a desired exact 
performance. No stochastic model can provide such rules as 
analytical models for adaptive algorithm behavior always rely 
on assumptions needed for mathematical tractability. Though 
the stationarity assumption is not satisfied in most practical 
systems, it is largely recognized that the models derived using 
them can still show tendencies of the algorithm behavior for 
reasonably small degrees of nonstationarity. To illustrate the 
validity of the model even for a nonstationary input signal, 
3 simulation scenarios were tested for a GSC-AEC system 
with M = 2, Ah = 1024, Abe = 16, Aaec = 1039, 
MAEC = Fbf = 3.1612 X 10“"^ and different degrees of 
nonstationarity rj n. Results are shown in Eig. [fi] Eor higher 
the degrees of nonstationarity, the behavior of the system 
diverges from the theoretical prediction. However, the system 
performance is still close enough to the MC simulation to 
jumpstart the design choices. 

D. Design Example 1 

Consider an acoustic echo cancellation system with a rever¬ 
beration time Tr( 60) = 60 ms (typical of a car cabin) and a 
background noise level of —20 dB. We assume the frequency 
response in the desired DOA (broadside) can be modeled 
by a linear phase, delayless all-pass filter with Ay = 16 
coefficients, i.e. / = [1, Oixatj--!]^ IMl- The design goals 
are convergence of J[n]dB < — 20 in less than 2 s (evaluated 
at n = 1.5A4 with (l62l i and (iTOl) '). 
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(b) rj = 0.1 


E 


MC simulation 
model prediction 


If more than one combination of (/iaeCi/ibf) is capable of 
reaching the desired cancellation at n = 1.5ii^4 then only the 
one with the lower J[oo] is considered. 



0.9). 


From Fig. |7] 2 candidate solutions were selected with 
J[cxd] Ri —21.5 dB and simulated using real speech signals 
(with pauses removed). Average results for an ensemble of 50 
runs are shown in Fig. |8] 


TABLE I 

Parameters eor Design Example 1 


M 

Naec 

P-AEC 


MBF 


J[cxd] 

2 

405 

9.7778E - 

04 

6.4603F; - 

04 

-21.54 dB 

4 

290 

9.1034F;- 

04 

1.6969F; - 

04 

-21.5 dB 


s -10 



.25 1- ^^^^- 

0 1 2 3 4 5 


" xio'^ 

(c) T) = 0.5 

Fig. 6. MC simulation for Model Verification 3 

We consider the design using the step-size matrix in (l44li . 
The frequency response model requires A^bf = Nj = IQ. The 
LEM plants have length > /sTr(60) = 480, and we thus 
set Ah = 500. The free parameters are then M and Waec- The 
choices of M and Waec affect the computational complexity 
per iteration and the convergence speed. 

Fig. m was produced evaluating (l62l i and (|70l for 
M = {1,2}, Naec from 290 to 515, and tr(A4i2mod) = 
(2/300,..., 2/3} to compute J[cx)]. If J[oo] > —20 dB this 
configuration is discarded. (l62] | and (TTOI i are evaluated for 
fJ-AEC = (0.01,... 0.99} X 2( J[oo] - Ji„in)/(J'[oo]tr(i2u.u. )). 



Fig. 8. MC simulation (speech signals) 


E. Design Example 2 

To illustrate the use of the derived optimal step-matrix de¬ 
rived in (ISTT i. consider an acoustic echo canceler with M = 2 
microphones for a large conference room with reverberation 
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time Tf>( 60) = 0.45 s. Assuming a sampling rate of 8 kHz, 
the LEM plant length is about Ah ~ 0.45.8 x 10^ = 3600 
coefficients. For this sampling rate, we also consider the 
frequency response in the desired DOA can be guaranteed 
with Nf = 16 constraints. In this simulation we consider only 
design choices with Aaec = Ah +Abf — 1. The far-end signal 
is modeled by a unity variance AR1(—0.9) random process and 
the noises in each microphone are assumed independent and 
modeled by Gaussian i.i.d. variables with variance 10“^. In 
this design we desire a steady-state MOP of —22 dB and a —10 
dB MOP after 1 second of convergence (n = 8000). To verify 
the feasibility of this design we used the proposed statistical 
model to predict the transient behavior using the optimal step- 
matrix derived in JsTl l. Results are shown in Fig.|9] We observe 



Fig. 9. Model prediction for Design Example 2 

that the whitening transformation is capable of achieving a 
MOP around —8 dB at n = 8000. As this is the design with the 
optimal convergence speed, it is not possible to achieve, at the 
same time, the desired steady-state MOP of —22 dB and less 
than —8 dB at n = 8000. Therefore it is not possible to design 
an BF-AEC system in the GSC form with Tfi(60) = 0.45 and 
these performance requirements. 

F. Design Example 3 

Consider an BF-AEC system designed to work in a room 
with reverberation time rR(60) = 100 ms. To guarantee a 
Public Switched Telephone Network (PSTN) quality signal, 
the sampling rate was chosen fs — 8000. To model the 
LEM plant impulse responses we have chosen Aj, = 1000 
coefficients. Assume there are M = 2 microphones available 
and a reasonable frequency response in the DOA is achieved 
with Abf = 16. Finally, an AFC length Aafc = .^h + 
^BF — 1 = 1015 coefficients was used. For this example, 
we assume an ideal double-talk detector is available. The 
initial DOA is assumed initially at a 7r/4 angle in relation 
to the broadside of the microphone array. During the first 
10® samples, the adaptation occurs in the absence of near-end 
speech with equal step-sizes paec = Mbf = 8.2147 x 10“®. 
Then, a double-talk period occurs. The unitary power near¬ 
end speech, modelled as an ARl(-0.9) process, arrives from 


the broadside of the microphone array. We assume double- 
talk control logic, constraint and block matrix correction act 
instantly. The adaptation of the AFC is frozen by the double- 
talk detector (p-aec = 0), and during the next 5.10® samples 
only the BF is adapted with /ibe = 8.2147 x 10“®. During this 
period, the convergence is significantly faster as the BF-AEC 
structure is not jointly-optimized and the effective adaptive 
filter length is reduced to MAbf — ^f- During the next 10® 
samples, the near-end speech is removed and, considering the 
BF state from the double-talk period is a good initial solution, 
step-sizes are set to /taec = 9.3243 x 10“® and pbf = 10”^. 
Finally, the LEM plant is subjected to an abrupt change in 
which a completely new H is used. In this configuration the 
new step-sizes are set to /taec = Mbf = 9x 10“® to accelerate 
convergence after an abrupt LEM plant change I?] . The model 
predictions and Monte Carlo Simulations (ensemble of 50 
runs) are compared on Fig. [TO] 



Fig. 10. Model prediction and Monte-Carlo simulation results Design 
Example 3 

IX. Conclusion 

This work presented a statistical analysis of a class of 
jointly optimized beamformer-assisted AFC. The analysis was 
performed for systems with the BF implemented in the GSC 
form and using the LMS algorithm. The analysis considered 
convergence control using a step size matrix to accommodate 
typical control logic implementations. We have shown that the 
joint optimization of the BF-AEC is equivalent to a LCMV 
problem. Thus, the derived analytical models can be used to 
predict the transient performance of general adaptive wideband 
beamformers. The stochastic model was determined for the 
transient and steady-state behaviors of the residual mean echo 
power for stationary Gaussian inputs. Convergence analysis 
lead to stability bounds for the step-size matrix. Design 
guidelines were derived from the analytical models. Monte 
Carlo simulations illustrated the accuracy of the theoretical 
models and the applicability of the proposed design guidelines. 
Finally, it was shown how a high convergence rate can be 
achieved using a quasi-Newton adaptation scheme in which 
the step-size matrix is designed to whiten the combined input 
vector. 
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